package Jama.examples;
import java.util.Date;
import Jama.EigenvalueDecomposition;
import Jama.LUDecomposition;
import Jama.Matrix;
import Jama.QRDecomposition;
/** Example of use of Matrix Class, featuring magic squares. **/
public class MagicSquareExample {
/** Generate magic square test matrix. **/
public static Matrix magic(int n) {
double[][] M = new double[n][n];
// Odd order
if ((n % 2) == 1) {
int a = (n + 1) / 2;
int b = (n + 1);
for (int j = 0; j < n; j++) {
for (int 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 (int j = 0; j < n; j++) {
for (int 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 {
int p = n / 2;
int k = (n - 2) / 4;
Matrix A = magic(p);
for (int j = 0; j < p; j++) {
for (int i = 0; i < p; i++) {
double 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 (int i = 0; i < p; i++) {
for (int j = 0; j < k; j++) {
double t = M[i][j];
M[i][j] = M[i + p][j];
M[i + p][j] = t;
}
for (int j = n - k + 1; j < n; j++) {
double t = M[i][j];
M[i][j] = M[i + p][j];
M[i + p][j] = t;
}
}
double 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);
}
/** Shorten spelling of print. **/
private static void print(String s) {
System.out.print(s);
}
/** Format double with Fw.d. **/
public static String fixedWidthDoubletoString(double x, int w, int d) {
java.text.DecimalFormat fmt = new java.text.DecimalFormat();
fmt.setMaximumFractionDigits(d);
fmt.setMinimumFractionDigits(d);
fmt.setGroupingUsed(false);
String s = fmt.format(x);
while (s.length() < w) {
s = " " + s;
}
return s;
}
/** Format integer with Iw. **/
public static String fixedWidthIntegertoString(int n, int w) {
String s = Integer.toString(n);
while (s.length() < w) {
s = " " + s;
}
return s;
}
public static void main(String argv[]) {
/*
| 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).
*/
print("\n Test of Matrix Class, using magic squares.\n");
print(" See MagicSquareExample.main() for an explanation.\n");
print(
"\n n trace max_eig rank cond lu_res qr_res\n\n");
Date start_time = new Date();
double eps = Math.pow(2.0, -52.0);
for (int n = 3; n <= 32; n++) {
print(fixedWidthIntegertoString(n, 7));
Matrix M = magic(n);
int t = (int) M.trace();
print(fixedWidthIntegertoString(t, 10));
EigenvalueDecomposition E = new EigenvalueDecomposition(
M.plus(M.transpose()).times(0.5));
double[] d = E.getRealEigenvalues();
print(fixedWidthDoubletoString(d[n - 1], 14, 3));
int r = M.rank();
print(fixedWidthIntegertoString(r, 7));
double c = M.cond();
print(
c < 1 / eps
? fixedWidthDoubletoString(c, 12, 3)
: " Inf");
LUDecomposition LU = new LUDecomposition(M);
Matrix L = LU.getL();
Matrix U = LU.getU();
int[] p = LU.getPivot();
Matrix R = L.times(U).minus(M.getMatrix(p, 0, n - 1));
double res = R.norm1() / (n * eps);
print(fixedWidthDoubletoString(res, 12, 3));
QRDecomposition QR = new QRDecomposition(M);
Matrix Q = QR.getQ();
R = QR.getR();
R = Q.times(R).minus(M);
res = R.norm1() / (n * eps);
print(fixedWidthDoubletoString(res, 12, 3));
print("\n");
}
Date stop_time = new Date();
double etime = (stop_time.getTime() - start_time.getTime()) / 1000.;
print(
"\nElapsed Time = " + fixedWidthDoubletoString(etime, 12, 3)
+ " seconds\n");
print("Adios\n");
}
}