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");
   }
}