Jama.examples.MagicSquareExample Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of jama Show documentation
Show all versions of jama Show documentation
JAMA is a basic linear algebra package for Java. It provides
user-level classes for constructing and manipulating real, dense matrices.
It is meant to provide sufficient functionality for routine problems,
packaged in a way that is natural and understandable to
non-experts.
package Jama.examples;
import Jama.*;
import java.util.Date;
/** 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");
}
}