|
/Users/lyon/j4p/src/math/linearAlgebra/examples/MagicSquareExample.java
|
1 package math.linearAlgebra.examples;
2
3 import math.linearAlgebra.EigenvalueDecomposition;
4 import math.linearAlgebra.LUDecomposition;
5 import math.linearAlgebra.Matrix;
6 import math.linearAlgebra.QRDecomposition;
7
8 import java.util.Date;
9
10 /** Example of use of Matrix Class, featuring magic squares. **/
11
12 public class MagicSquareExample {
13
14 /** Generate magic square test matrix. **/
15
16 public static Matrix magic(int n) {
17
18 double[][] M = new double[n][n];
19
20 // Odd order
21
22 if ((n % 2) == 1) {
23 int a = (n + 1) / 2;
24 int b = (n + 1);
25 for (int j = 0; j < n; j++) {
26 for (int i = 0; i < n; i++) {
27 M[i][j] = n * ((i + j + a) % n) + ((i + 2 * j + b) % n) + 1;
28 }
29 }
30
31 // Doubly Even Order
32
33 } else if ((n % 4) == 0) {
34 for (int j = 0; j < n; j++) {
35 for (int i = 0; i < n; i++) {
36 if (((i + 1) / 2) % 2 == ((j + 1) / 2) % 2) {
37 M[i][j] = n * n - n * i - j;
38 } else {
39 M[i][j] = n * i + j + 1;
40 }
41 }
42 }
43
44 // Singly Even Order
45
46 } else {
47 int p = n / 2;
48 int k = (n - 2) / 4;
49 Matrix A = magic(p);
50 for (int j = 0; j < p; j++) {
51 for (int i = 0; i < p; i++) {
52 double aij = A.get(i, j);
53 M[i][j] = aij;
54 M[i][j + p] = aij + 2 * p * p;
55 M[i + p][j] = aij + 3 * p * p;
56 M[i + p][j + p] = aij + p * p;
57 }
58 }
59 for (int i = 0; i < p; i++) {
60 for (int j = 0; j < k; j++) {
61 double t = M[i][j];
62 M[i][j] = M[i + p][j];
63 M[i + p][j] = t;
64 }
65 for (int j = n - k + 1; j < n; j++) {
66 double t = M[i][j];
67 M[i][j] = M[i + p][j];
68 M[i + p][j] = t;
69 }
70 }
71 double t = M[k][0];
72 M[k][0] = M[k + p][0];
73 M[k + p][0] = t;
74 t = M[k][k];
75 M[k][k] = M[k + p][k];
76 M[k + p][k] = t;
77 }
78 return new Matrix(M);
79 }
80
81 /** Shorten spelling of print. **/
82
83 private static void print(String s) {
84 System.out.print(s);
85 }
86
87 /** Format double with Fw.d. **/
88
89 public static String fixedWidthDoubletoString(double x, int w, int d) {
90 java.text.DecimalFormat fmt = new java.text.DecimalFormat();
91 fmt.setMaximumFractionDigits(d);
92 fmt.setMinimumFractionDigits(d);
93 fmt.setGroupingUsed(false);
94 String s = fmt.format(x);
95 while (s.length() < w) {
96 s = " " + s;
97 }
98 return s;
99 }
100
101 /** Format integer with Iw. **/
102
103 public static String fixedWidthIntegertoString(int n, int w) {
104 String s = Integer.toString(n);
105 while (s.length() < w) {
106 s = " " + s;
107 }
108 return s;
109 }
110
111
112 public static void main(String argv[]) {
113
114 /*
115 | Tests LU, QR, SVD and symmetric Eig decompositions.
116 |
117 | n = order of magic square.
118 | trace = diagonal sum, should be the magic sum, (n^3 + n)/2.
119 | max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
120 | rank = linear algebraic rank,
121 | should equal n if n is odd, be less than n if n is even.
122 | cond = L_2 condition number, ratio of singular values.
123 | lu_res = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
124 | qr_res = test of QR factorization, norm1(Q*R-A)/(n*eps).
125 */
126
127 print("\n Test of Matrix Class, using magic squares.\n");
128 print(" See MagicSquareExample.main() for an explanation.\n");
129 print("\n n trace max_eig rank cond lu_res qr_res\n\n");
130
131 Date start_time = new Date();
132 double eps = Math.pow(2.0, -52.0);
133 for (int n = 3; n <= 32; n++) {
134 print(fixedWidthIntegertoString(n, 7));
135
136 Matrix M = magic(n);
137
138 int t = (int) M.trace();
139 print(fixedWidthIntegertoString(t, 10));
140
141 EigenvalueDecomposition E =
142 new EigenvalueDecomposition(M.plus(M.transpose()).times(0.5));
143 double[] d = E.getRealEigenvalues();
144 print(fixedWidthDoubletoString(d[n - 1], 14, 3));
145
146 int r = M.rank();
147 print(fixedWidthIntegertoString(r, 7));
148
149 double c = M.cond();
150 print(c < 1 / eps ? fixedWidthDoubletoString(c, 12, 3) :
151 " Inf");
152
153 LUDecomposition LU = new LUDecomposition(M);
154 Matrix L = LU.getL();
155 Matrix U = LU.getU();
156 int[] p = LU.getPivot();
157 Matrix R = L.times(U).minus(M.getMatrix(p, 0, n - 1));
158 double res = R.norm1() / (n * eps);
159 print(fixedWidthDoubletoString(res, 12, 3));
160
161 QRDecomposition QR = new QRDecomposition(M);
162 Matrix Q = QR.getQ();
163 R = QR.getR();
164 R = Q.times(R).minus(M);
165 res = R.norm1() / (n * eps);
166 print(fixedWidthDoubletoString(res, 12, 3));
167
168 print("\n");
169 }
170 Date stop_time = new Date();
171 double etime = (stop_time.getTime() - start_time.getTime()) / 1000.;
172 print("\nElapsed Time = " +
173 fixedWidthDoubletoString(etime, 12, 3) + " seconds\n");
174 print("Adios\n");
175 }
176 }
177