<?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();
?>
