<?php /** * @package JAMA */ require_once "../Matrix.php"; /** * Example of use of Matrix Class, featuring magic squares. */ class MagicSquareExample { /** * Generate magic square test matrix. * @param int n dimension of matrix */ function magic($n) { // Odd order if (($n % 2) == 1) { $a = ($n+1)/2; $b = ($n+1); for ($j = 0; $j < $n; ++$j) for ($i = 0; $i < $n; ++$i) $M[$i][$j] = $n*(($i+$j+$a) % $n) + (($i+2*$j+$b) % $n) + 1; // Doubly Even Order } else if (($n % 4) == 0) { for ($j = 0; $j < $n; ++$j) { for ($i = 0; $i < $n; ++$i) { if ((($i+1)/2)%2 == (($j+1)/2)%2) $M[$i][$j] = $n*$n-$n*$i-$j; else $M[$i][$j] = $n*$i+$j+1; } } // Singly Even Order } else { $p = $n/2; $k = ($n-2)/4; $A = $this->magic($p); $M = array(); for ($j = 0; $j < $p; ++$j) { for ($i = 0; $i < $p; ++$i) { $aij = $A->get($i,$j); $M[$i][$j] = $aij; $M[$i][$j+$p] = $aij + 2*$p*$p; $M[$i+$p][$j] = $aij + 3*$p*$p; $M[$i+$p][$j+$p] = $aij + $p*$p; } } for ($i = 0; $i < $p; ++$i) { for ($j = 0; $j < $k; ++$j) { $t = $M[$i][$j]; $M[$i][$j] = $M[$i+$p][$j]; $M[$i+$p][$j] = $t; } for ($j = $n-$k+1; $j < $n; ++$j) { $t = $M[$i][$j]; $M[$i][$j] = $M[$i+$p][$j]; $M[$i+$p][$j] = $t; } } $t = $M[$k][0]; $M[$k][0] = $M[$k+$p][0]; $M[$k+$p][0] = $t; $t = $M[$k][$k]; $M[$k][$k] = $M[$k+$p][$k]; $M[$k+$p][$k] = $t; } return new Matrix($M); } /** * Simple function to replicate PHP 5 behaviour */ function microtime_float() { list($usec, $sec) = explode(" ", microtime()); return ((float)$usec + (float)$sec); } /** * Tests LU, QR, SVD and symmetric Eig decompositions. * * n = order of magic square. * trace = diagonal sum, should be the magic sum, (n^3 + n)/2. * max_eig = maximum eigenvalue of (A + A')/2, should equal trace. * rank = linear algebraic rank, should equal n if n is odd, * be less than n if n is even. * cond = L_2 condition number, ratio of singular values. * lu_res = test of LU factorization, norm1(L*U-A(p,:))/(n*eps). * qr_res = test of QR factorization, norm1(Q*R-A)/(n*eps). */ function main() { ?> <p>Test of Matrix Class, using magic squares.</p> <p>See MagicSquareExample.main() for an explanation.</p> <table border='1' cellspacing='0' cellpadding='4'> <tr> <th>n</th> <th>trace</th> <th>max_eig</th> <th>rank</th> <th>cond</th> <th>lu_res</th> <th>qr_res</th> </tr> <?php $start_time = $this->microtime_float(); $eps = pow(2.0,-52.0); for ($n = 3; $n <= 6; ++$n) { echo "<tr>"; echo "<td align='right'>$n</td>"; $M = $this->magic($n); $t = (int) $M->trace(); echo "<td align='right'>$t</td>"; $O = $M->plus($M->transpose()); $E = new EigenvalueDecomposition($O->times(0.5)); $d = $E->getRealEigenvalues(); echo "<td align='right'>".$d[$n-1]."</td>"; $r = $M->rank(); echo "<td align='right'>".$r."</td>"; $c = $M->cond(); if ($c < 1/$eps) echo "<td align='right'>".sprintf("%.3f",$c)."</td>"; else echo "<td align='right'>Inf</td>"; $LU = new LUDecomposition($M); $L = $LU->getL(); $U = $LU->getU(); $p = $LU->getPivot(); // Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1)); $S = $L->times($U); $R = $S->minus($M->getMatrix($p,0,$n-1)); $res = $R->norm1()/($n*$eps); echo "<td align='right'>".sprintf("%.3f",$res)."</td>"; $QR = new QRDecomposition($M); $Q = $QR->getQ(); $R = $QR->getR(); $S = $Q->times($R); $R = $S->minus($M); $res = $R->norm1()/($n*$eps); echo "<td align='right'>".sprintf("%.3f",$res)."</td>"; echo "</tr>"; } echo "<table>"; echo "<br />"; $stop_time = $this->microtime_float(); $etime = $stop_time - $start_time; echo "<p>Elapsed time is ". sprintf("%.4f",$etime) ." seconds.</p>"; } } $magic = new MagicSquareExample(); $magic->main(); ?>