File Source for MagicSquareExample.php

JAMA

# Source for file MagicSquareExample.php

Documentation is available at MagicSquareExample.php

1. <?php
2. /**
3. * @package JAMA
4. */
6. require_once "../Matrix.php";
8. /**
9. * Example of use of Matrix Class, featuring magic squares.
10. */
11. class MagicSquareExample {
13. /**
14. * Generate magic square test matrix.
15. * @param int n dimension of matrix
16. */
17. function magic($n) {
19. // Odd order
21. if (($n % 2) == 1) {
22. $a = ($n+1)/2;
23. $b = ($n+1);
24. for ($j = 0; $j < $n; ++$j)
25. for ($i = 0; $i < $n; ++$i)
26. $M[$i][$j] = $n*(($i+$j+$a) % $n) + (($i+2*$j+$b) % $n) + 1;
28. // Doubly Even Order
30. } else if (($n % 4) == 0) {
31. for ($j = 0; $j < $n; ++$j) {
32. for ($i = 0; $i < $n; ++$i) {
33. if ((($i+1)/2)%2 == (($j+1)/2)%2)
34. $M[$i][$j] = $n*$n-$n*$i-$j;
35. else
36. $M[$i][$j] = $n*$i+$j+1;
37. }
38. }
40. // Singly Even Order
42. } else {
44. $p = $n/2;
45. $k = ($n-2)/4;
46. $A = $this-> magic($p);
47. $M = array();
48. for ($j = 0; $j < $p; ++$j) {
49. for ($i = 0; $i < $p; ++$i) {
50. $aij = $A->get($i,$j);
51. $M[$i][$j] = $aij;
52. $M[$i][$j+$p] = $aij + 2*$p*$p;
53. $M[$i+$p][$j] = $aij + 3*$p*$p;
54. $M[$i+$p][$j+$p] = $aij + $p*$p;
55. }
56. }
58. for ($i = 0; $i < $p; ++$i) {
59. for ($j = 0; $j < $k; ++$j) {
60. $t = $M[$i][$j];
61. $M[$i][$j] = $M[$i+$p][$j];
62. $M[$i+$p][$j] = $t;
63. }
64. for ($j = $n-$k+1; $j < $n; ++$j) {
65. $t = $M[$i][$j];
66. $M[$i][$j] = $M[$i+$p][$j];
67. $M[$i+$p][$j] = $t;
68. }
69. }
71. $t = $M[$k][0]; $M[$k][0] = $M[$k+$p][0]; $M[$k+$p][0] = $t;
72. $t = $M[$k][$k]; $M[$k][$k] = $M[$k+$p][$k]; $M[$k+$p][$k] = $t;
74. }
76. return new Matrix($M);
78. }
80. /**
81. * Simple function to replicate PHP 5 behaviour
82. */
83. function microtime_float() {
84. list($usec, $sec) = explode(" ", microtime());
85. return ((float)$usec + (float)$sec);
86. }
88. /**
89. * Tests LU, QR, SVD and symmetric Eig decompositions.
90. *
91. * n = order of magic square.
92. * trace = diagonal sum, should be the magic sum, (n^3 + n)/2.
93. * max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
94. * rank = linear algebraic rank, should equal n if n is odd,
95. * be less than n if n is even.
96. * cond = L_2 condition number, ratio of singular values.
97. * lu_res = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
98. * qr_res = test of QR factorization, norm1(Q*R-A)/(n*eps).
99. */
100. function main() {
101. ?>
102. <p>Test of Matrix Class, using magic squares.</p>
103. <p>See MagicSquareExample.main() for an explanation.</p>
104. <table border='1' cellspacing='0' cellpadding='4'>
105. <tr>
106. <th>n</th>
107. <th>trace</th>
108. <th>max_eig</th>
109. <th>rank</th>
110. <th>cond</th>
111. <th>lu_res</th>
112. <th>qr_res</th>
113. </tr>
114. <?php
115. $start_time = $this-> microtime_float();
116. $eps = pow(2.0,-52.0);
117. for ($n = 3; $n <= 6; ++$n) {
118. echo "<tr>";
120. echo "<td align='right'>$n</td>";
122. $M = $this-> magic($n);
123. $t = (int) $M->trace();
125. echo "<td align='right'>$t</td>";
127. $O = $M->plus($M->transpose());
128. $E = new EigenvalueDecomposition($O->times(0.5));
129. $d = $E->getRealEigenvalues();
131. echo "<td align='right'>".$d[$n-1]."</td>";
133. $r = $M->rank();
135. echo "<td align='right'>".$r."</td>";
137. $c = $M->cond();
139. if ($c < 1/$eps)
140. echo "<td align='right'>".sprintf("%.3f",$c)."</td>";
141. else
142. echo "<td align='right'>Inf</td>";
144. $LU = new LUDecomposition($M);
145. $L = $LU->getL();
146. $U = $LU->getU();
147. $p = $LU->getPivot();
148. // Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1));
149. $S = $L->times($U);
150. $R = $S->minus($M->getMatrix($p,0,$n-1));
151. $res = $R->norm1()/($n*$eps);
153. echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
155. $QR = new QRDecomposition($M);
156. $Q = $QR->getQ();
157. $R = $QR->getR();
158. $S = $Q->times($R);
159. $R = $S->minus($M);
160. $res = $R->norm1()/($n*$eps);
162. echo "<td align='right'>".sprintf("%.3f",$res)."</td>";
164. echo "</tr>";
166. }
167. echo "<table>";
168. echo "<br />";
170. $stop_time = $this-> microtime_float();
171. $etime = $stop_time - $start_time;
173. echo "<p>Elapsed time is ". sprintf("%.4f",$etime) ." seconds.</p>";
175. }
177. }
179. $magic = new MagicSquareExample();
180. $magic-> main();
182. ?>

Documentation generated on Sun, 27 Feb 2011 16:32:34 -0800 by phpDocumentor 1.4.3