share/contrib/diag_test.mac

   1 /*
   2
   3   Simple testsuite for diag.mac
   4
   5   Run this file with
   6     batch ("diag_test", 'test)
   7 */
   8 kill (all);
   9 done$
  10
  11 (load ("diag"), 0);
  12 0$
  13
  14 diag ([1,2,3]);
  15 matrix ([1,0,0], [0,2,0], [0,0,3])$
  16
  17 diag ([matrix ([1]), matrix ([2]), matrix ([3])]);
  18 matrix ([1,0,0], [0,2,0], [0,0,3])$
  19
  20 (a1: matrix ([1,2,3], [0,4,5], [0,0,6]),
  21  a2: matrix([1,1], [1,0]),
  22  0);
  23 0$
  24
  25 diag ([a1, x, a2]);
  26 matrix ([1,2,3,0,0,0],
  27         [0,4,5,0,0,0],
  28         [0,0,6,0,0,0],
  29         [0,0,0,x,0,0],
  30         [0,0,0,0,1,1],
  31         [0,0,0,0,1,0])$
  32
  33 kill (a1, a2);
  34 done$
  35
  36 JF (0, 1);
  37 matrix ([0])$
  38
  39 JF (x, 3);
  40 matrix ([x, 1, 0],
  41         [0, x, 1],
  42         [0, 0, x])$
  43
  44 dispJordan ([[0, 1], [2, 2], [3, 2, 1]]);
  45 ''(diag ([JF (0, 1), JF (2, 2), JF (3, 2), JF (3, 1)]))$
  46
  47 /*
  48   For testing jordan(), use JCF examples built by dispJordan and
  49   conjugate by invertible matrices.
  50 */
  51 (make_conjugate (j_info) := block (
  52     [n: lsum (x, x, map (lambda ([lst], lsum (x, x, rest (lst))), j_info)),
  53      T],
  54     T: genmatrix (lambda([i,j], if is (i = n+1-j) then 0 else 1), n, n),
  55     T^^(-1) . dispJordan (j_info) . T),
  56  test_jordan (j_info) := is (jordan (make_conjugate (j_info)) = j_info),
  57  test_modematrix (j_info) := block (
  58    [conj: make_conjugate (j_info), T],
  59    T: ModeMatrix (conj, j_info),
  60    is (T^^(-1) . conj . T = dispJordan (j_info))),
  61  0);
  62 0$
  63
  64 jordan (matrix ([1]));
  65 [[1, 1]]$
  66
  67 test_jordan ([[%i, 2]]);
  68 true$
  69
  70 test_jordan ([[1, 1], [2, 2, 1]]);
  71 true$
  72
  73 minimalPoly ([[%i, 1]]);
  74 x - %i$
  75
  76 minimalPoly ([[1, 2, 1]]);
  77 (x - 1)^2$
  78
  79 minimalPoly ([[1, 1], [2, 1]]);
  80 (x - 1)*(x - 2)$
  81
  82 diag_sorted_list_histogram ([]);
  83 []$
  84
  85 diag_sorted_list_histogram ([1]);
  86 [[1, 1]]$
  87
  88 diag_sorted_list_histogram ([3, 3, 2, 2, 1]);
  89 [[3, 2], [2, 2], [1, 1]]$
  90
  91 /*
  92   Testing diag_kernel_element, diag_general_jordan_chain and
  93   diag_find_li_chain is probably a bit difficult. Let's punt on that
  94   and hope they're well exercised by ModeMatrix.
  95 */
  96 test_modematrix ([[%i, 2]]);
  97 true$
  98
  99 test_modematrix ([[1, 1], [2, 2, 1]]);
 100 true$
 101
 102 diag_taylor_coefficients (exp(x), x, 3);
 103 [exp(x), exp(x), exp(x)/2, exp(x)/6]$
 104
 105 diag_taylor_expand_block (diag_taylor_coefficients (exp(x), x, 3),
 106                           x, 1, 3);
 107 matrix ([exp(1), exp(1), exp(1)/2],
 108         [0,      exp(1), exp(1)],
 109         [0,      0,      exp(1)])$
 110
 111 mat_function (sin, ident(3)*t);
 112 ''(ident(3)*sin(t))$
 113
 114 /*
 115   Bug 2542
 116
 117   The "correct" answer was calculated numerically by Octave. There's
 118   probably a better way of checking this!
 119 */
 120 block ([a: matrix([6,6,-3,2],[-4,-4,2,0],[8,7,-4,4],[1,0,-1,-2]),
 121         num_ans:
 122           matrix ([ 2.9848e+00,  2.3771e+00, -1.2382e+00,  6.1912e-01],
 123                   [-1.2382e+00, -7.2991e-01,  6.1912e-01,  1.3645e-16],
 124                   [ 3.7147e+00,  3.2064e+00, -1.3490e+00,  1.2382e+00],
 125                   [-1.1079e-01, -2.6557e-01, -9.9383e-02, -1.1079e-01]),
 126         err],
 127   err: abs (float (rectform (mat_function (exp, a))) - num_ans),
 128   is (lmax (map (lmax, args (err))) < 1e-4));
 129 true$