I stumbled onto a library for manipulating and analyzing matrices. I’m pretty sure there are several options for doing this in Processing, but this option here works so I thought I would share the example.
First, save a sketch and then add a “code” directory. Download the Jama(java matrix library) jar file from this link: https://math.nist.gov/javanumerics/jama/
Then put this code into the saved sketch and save the downloaded jar file into the code directory of the sketch. If you download the zip file version, you get the source code and the documentation for the library.:
import Jama.*;
//import controlP5.*;
//import Jama.Matrix;
import java.util.Date;
//class MagicSquareExample {
/** Generate magic square test matrix. **/
public 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);
}
public 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;
}
public String fixedWidthIntegertoString (int n, int w) {
String s = Integer.toString(n);
while (s.length() < w) {
s = " " + s;
}
return s;
}
//ControlP5 cp5;
//ControlP5 controlP5;
//Textarea myTextarea;
void setup() {
// fullscreen();
size(1275, 750, OPENGL);
/* myTextarea = cp5.addTextarea("txt")
.setPosition(450,10)
.setSize(800,700)
.setFont(createFont("arial",14))
.setLineHeight(14)
.setColor(color(255))
.setColorBackground(color(255,0))
.setColorForeground(color(255,100)); */
}
void draw() {
// public void main (String argv[]) {
//noLoop();
/*
| 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");
}