share/linearalgebra/rtest_lu.mac

   1 (kill(values),0);
   2 0$
   3
   4 (ratprint : false, float2bf : true,0);
   5 0$
   6
   7 (lu_test(m) := block([mp],
   8   mp : get_lu_factors(lu_factor(m,generalring)),
   9   zeromatrixp(m - mp[1] . mp[2] . mp[3])),0);
  10 0$
  11
  12 lu_test(matrix([1]));
  13 true$
  14
  15 lu_test(matrix([0,1],[1,0]));
  16 true$
  17
  18 lu_test(matrix([0,a],[b,0]));
  19 true$
  20
  21 lu_test(matrix([a,b],[c,d]));
  22 true$
  23
  24 lu_test(matrix([a,b],[b,d]));
  25 true$
  26
  27 lu_test(matrix([1/2,b],[b,2/3]));
  28 true$
  29
  30 lu_test(matrix([1/2 + %i, %i - 1],[9 - %i, 19 + %i]));
  31 true$
  32
  33 lu_test(matrix([0,1,0],[0,0,1],[1,2,4]));
  34 true$
  35
  36 lu_test(matrix([0,a,0],[b,0,0],[0,0,c]));
  37 true$
  38
  39 lu_test(matrix([a,b,c],[d,e,f],[g,h,i]));
  40 true$
  41
  42 lu_test(matrix([0,0,1],[1,0,0],[0,1,0]));
  43 true$
  44
  45 lu_test(matrix([0,0,1],[1,0,0],[0,1,x]));
  46 true$
  47
  48 lu_test(hilbert_matrix(3));
  49 true$
  50
  51 lu_test(hilbert_matrix(4));
  52 true$
  53
  54 lu_test(hilbert_matrix(5));
  55 true$
  56
  57 lu_test(vandermonde_matrix([1,2,%i]));
  58 true$
  59
  60 lu_test(vandermonde_matrix([1 + %i,1 - %i, 8, 32]));
  61 true$
  62
  63 (lu_test(m, tol) := block([mp],
  64   mp : get_lu_factors(lu_factor(m,bigfloatfield)),
  65   every(lambda([s], is(cabs(s) < tol)), m - mp[1] . mp[2] . mp[3])),0);
  66 0$
  67
  68 (fpprec : 20,0);
  69 0$
  70
  71 lu_test(vandermonde_matrix([1,1/2,1/3,1/4]), 1.0e-18);
  72 true$
  73
  74 (fpprec : 40,0);
  75 0$
  76
  77 lu_test(vandermonde_matrix([1,1/2,1/3,1/4]), 1.0e-38);
  78 true$
  79
  80 lu_test(hilbert_matrix(13), 1.0e-38);
  81 true$
  82
  83 lu_test(hilbert_matrix(4) + %i * vandermonde_matrix([1,2,3,4]), 1.0e-15);
  84 true$
  85
  86 (lu_test(m, tol) := block([mp],
  87   mp : get_lu_factors(lu_factor(m,complexfield)),
  88   every(lambda([s], is(cabs(s) < tol)),m - mp[1] . mp[2] . mp[3])),0);
  89 0$
  90
  91 lu_test(hilbert_matrix(7) + %i * vandermonde_matrix([1,2,3,4,5,6,7]), 1.0e-10);
  92 true$
  93
  94 (m : matrix([1,2,3],[4,5,6],[7,8,10]),0);
  95 0$
  96
  97 (b : m . matrix([x],[y],[z]),0);
  98 0$
  99
 100 lu_backsub(lu_factor(m,crering), b);
 101 ''(rat(matrix([x],[y],[z])))$
 102
 103 (m : matrix([0,2,3],[0,5,6],[7,8,10]),0);
 104 0$
 105
 106 (b : m . matrix([x],[y],[z]),0);
 107 0$
 108
 109 lu_backsub(lu_factor(m, crering), b);
 110 ''(rat(matrix([x],[y],[z])))$
 111
 112 (m : matrix([0,2,3],[0,5,6],[7,8,10]),0);
 113 0$
 114
 115 (b : m . matrix([x],[y],[z]),0);
 116 0$
 117
 118 lu_backsub(lu_factor(m, crering), b);
 119 ''(rat(matrix([x],[y],[z])))$
 120
 121
 122 (m : matrix([0,1,2],[1,0,1],[0,0,8]),0);
 123 0$
 124
 125 (b : m . matrix([x],[y],[z]),0);
 126 0$
 127
 128 lu_backsub(lu_factor(m, crering), b);
 129 ''(rat(matrix([x],[y],[z])))$
 130
 131 (m : matrix([matrix([0,1],[1,0]), matrix([1,0],[%i,0])],
 132             [matrix([1,0],[0,1]), matrix([1,1],[-1,1])]),0);
 133 0$
 134
 135 (m1 : get_lu_factors(lu_factor(m, noncommutingring)),0);
 136 0$
 137
 138 (matrix_element_mult : ".",0);
 139 0$
 140
 141 zeromatrixp(m1[1] . m1[2] . m1[3] - m);
 142 true$
 143
 144 (m : matrix([matrix([matrix([8,4],[4,13]),matrix([2,4],[7,9])],[matrix([2,7],[4,9]),matrix([10,6],[6,15])]),matrix([matrix([4,4],[12,5]),matrix([0,6],[6,8])],[matrix([8,3],[13,3]),matrix([6,5],[5,6])])],[matrix([matrix([4,12],[4,5]),matrix([8,13],[3,3])],[matrix([0,6],[6,8]),matrix([6,5],[5,6])]),matrix([matrix([18,5],[5,4]),matrix([5,9],[2,6])],[matrix([5,2],[9,6]),matrix([6,3],[3,11])])]),0);
 145 0$
 146
 147 (m1 : get_lu_factors(lu_factor(m,'noncommutingring)),0);
 148 0$
 149
 150 zeromatrixp(m1[1] . m1[2] . m1[3]-m);
 151 true$
 152
 153 (matrix_element_mult : "*",0);
 154 0$
 155
 156 errcatch(lu_factor(matrix()));
 157 []$
 158
 159 get_lu_factors(lu_factor(matrix([0])));
 160 [matrix([1]),matrix([1]),matrix([0])]$
 161
 162 ratdisrep(get_lu_factors(lu_factor(matrix([0]),'crering)));
 163 [matrix([1]),matrix([1]),matrix([0])]$
 164
 165 get_lu_factors(lu_factor(matrix([0]), 'floatfield));
 166 [matrix([1]),matrix([1]),matrix([0.0])]$
 167
 168 get_lu_factors(lu_factor(matrix([0]), 'bigfloatfield));
 169 [matrix([1]),matrix([1]),matrix([0.0b0])]$
 170
 171 get_lu_factors(lu_factor(matrix([0,1],[0,1])));
 172 [matrix([0,1],[1,0]),matrix([1,0],[0,1]),matrix([0,1],[0,1])]$
 173
 174 get_lu_factors(lu_factor(matrix([5,0,0],[0,0,0],[0,0,1])));
 175 [matrix([1,0,0],[0,0,1],[0,1,0]),matrix([1,0,0],[0,1,0],[0,0,1]),
 176 matrix([5,0,0],[0,0,1],[0,0,0])]$
 177
 178 get_lu_factors(lu_factor(matrix([5,0,0],[0,0,0],[0,0,0])));
 179 [matrix([1,0,0],[0,0,1],[0,1,0]),
 180 matrix([1,0,0],[0,1,0],[0,0,1]),
 181 matrix([5,0,0],[0,0,0],[0,0,0])]$
 182
 183 (lower_triangular_p(mat) := block([n : matrix_size(mat), OK : true],
 184     for r :1 thru second(n) do (
 185       for c : r thru first(n) do (
 186             OK : OK and (if r=c then is (mat[r,c] = 1) else is(mat[r,c] = 0)))),
 187     OK),0);
 188 0$
 189
 190 (upper_triangular_p(mat) := block([n : matrix_size(mat), OK : true],
 191     for r :1 thru second(n) do (
 192       for c : 1 thru r-1 do (
 193             OK : OK and is(mat[r,c] = 0))),
 194     OK), 0);
 195 0$
 196
 197 (perm_mat_p(mat) := zeromatrixp(mat.transpose(mat) - identfor(mat)),0);
 198 0$
 199
 200 /* In this file, there are multiple definitions of lu_test." */
 201
 202
 203 (lu_test(m) := block([mp],
 204   mp : get_lu_factors(lu_factor(m, 'generalring)),
 205   lower_triangular_p(second(mp)) and
 206   perm_mat_p(first(mp)) and
 207   upper_triangular_p(third(mp)) and
 208   zeromatrixp(ratsimp(m - first(mp).second(mp).third(mp)))), 0);
 209 0$
 210
 211 /* see mailing list, "lu_factor fails on trivial diagonal matrix," 5 Jan 2021. */
 212 lu_test(matrix([-1,0,0,0],[0,-3,0,0],[0,0,0,0],[0,0,0,-3]));
 213 true$
 214
 215 lu_test(matrix([0,0,0,0],[0,-3,0,0],[0,0,0,0],[0,0,0,-3]));
 216 true$
 217
 218 lu_test(matrix([0,0,0,0],[0,0,0,0],[0,0,-3,0],[0,0,0,0]));
 219 true$
 220
 221 lu_test(matrix([0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,-3]));
 222 true$
 223
 224 lu_test(matrix([0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,-1]));
 225 true$
 226
 227 every(lambda([x],x), map('lu_test, flatten(outermap(lambda([a,b,c,d], matrix([a,b],[c,d])),
 228     [0,1],[0,1],[0,1],[0,1]))));
 229 true$
 230
 231 (remvalue(m,b,m1),0);
 232 0$
 233
 234 every(lambda([x],x), map(lu_test, flatten(outermap(lambda([x1,x2,x3,x4], matrix([x1,x2],[x3,x4])),
 235    [a,b],[a,b],[a,b],[a,b]))));
 236 true$
 237
 238 /* diagonal matrices */
 239 lu_test(matrix([46,0,0],[0,107,0],[0,0,28]));
 240 true$
 241
 242 lu_test(matrix([46,0,0],[0,107,0],[0,0,0]));
 243 true$
 244
 245 lu_test(matrix([0,0,0],[0,46,0],[0,0,107]));
 246 true$
 247
 248 lu_test(matrix([107,0,0],[0,0,0],[0,0,46]));
 249 true$
 250
 251 lu_test(matrix([0,0,0],[0,0,0],[0,0,46]));
 252 true$
 253
 254 lu_test(matrix([0,0,0],[0,107,0],[0,0,0]));
 255 true$
 256
 257 lu_test(matrix([107,0,0],[0,0,0],[0,0,0]));
 258 true$
 259
 260 lu_test(matrix([0,0,0],[0,0,0],[0,0,107]));
 261 true$
 262
 263 lu_test(matrix([0,0,0],[0,0,0],[0,0,0]));
 264 true$
 265
 266 lu_test(matrix([107,0,0],[0,46,0],[0,0,107]));
 267 true$
 268
 269 (mat : matrix([matrix([0,0],[0,1]), matrix([1,2],[2,1])], [matrix([0,1],[1,0]), matrix([1,2],[2,1])]),0);
 270 0$
 271
 272 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 273   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 274 true$
 275
 276 (mat : matrix([matrix([0,0],[0,0]), matrix([107,2],[2,107])], [matrix([0,0],[0,0]), matrix([1,2],[2,1])]),0);
 277 0$
 278
 279 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 280   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 281 true$
 282
 283 (mat : matrix([matrix([0,1],[1,0]), matrix([107,2],[2,107])], [matrix([0,0],[0,0]), matrix([1,2],[2,1])]),0);
 284 0$
 285
 286 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 287   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 288 true$
 289
 290 (mat : matrix([matrix([0,1],[1,0]), matrix([107,2],[2,107])], [matrix([0,1],[1,0]), matrix([1,46],[46,1])]),0);
 291 0$
 292
 293 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 294   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 295 true$
 296
 297 (mat : matrix([matrix([0,0],[1,0]), matrix([1,0],[0,0])], [matrix([0,1],[1,0]), matrix([1,46],[46,1])]),0);
 298 0$
 299
 300 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 301   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 302 true$
 303
 304 (mat : matrix([matrix([0,0],[0,0]), matrix([0,0],[0,0])], [matrix([0,0],[0,0]), matrix([0,0],[0,0])]),0);
 305 0$
 306
 307 block([matrix_element_mult : ".", mm : get_lu_factors(lu_factor(mat, 'noncommutingring))],
 308   zeromatrixp(mat - first(mm) . second(mm) . third(mm)));
 309 true$
 310
 311 (reset(ratprint, float2bf, matrix_element_mult, matrix_element_transpose),0);
 312 0$
 313
 314 (remfunction(lu_test, lower_triangular_p,upper_triangular_p,perm_mat_p),0);
 315 0$
 316
 317 (kill(values),0);
 318 0$
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329