share/linearalgebra/rtest_cholesky.mac

   1 (checkcholesky(m,r) := (
   2   l : cholesky(m,r),
   3   l : m - l . ctranspose(l),
   4   if blockmatrixp(l) then l : mat_fullunblocker(l),
   5   every(lambda([x], is(cabs(x) < 1.0e-10)),l)),0);
   6 0$
   7
   8 /*-----------*/
   9 (m : matrix([0]),0);
  10 0$
  11
  12 (ok : true,0);
  13 0$
  14
  15 (for r in [floatfield, complexfield, generalring, rationalfield, crering, noncommutingring] do (
  16   ok : ok and checkcholesky(m,r)),ok);
  17
  18 true$
  19
  20
  21 /*-----------*/
  22 (m : matrix([1]),0);
  23 0$
  24
  25 (ok : true,0);
  26 0$
  27
  28 (for r in [floatfield, complexfield, generalring, rationalfield, crering, noncommutingring] do (
  29   ok : ok and checkcholesky(m,r)),ok);
  30
  31 true$
  32
  33 /*-----------*/
  34 (m : matrix([1,0],[7,-8]),0);
  35 0$
  36
  37 (m : m . ctranspose(m),0);
  38 0$
  39
  40 (ok : true,0);
  41 0$
  42
  43 (for r in [floatfield, complexfield, generalring, crering, noncommutingring] do (
  44   ok : ok and checkcholesky(m,r)),ok);
  45 true$
  46
  47 /*-----------*/
  48 (m : matrix([1,0],[7,-%i]),0);
  49 0$
  50
  51 (m : m . ctranspose(m),0);
  52 0$
  53
  54 (ok : true,0);
  55 0$
  56
  57 (for r in [complexfield, generalring] do (
  58   ok : ok and checkcholesky(m,r)),ok);
  59 true$
  60
  61
  62 /*-----------*/
  63 (m : matrix([1,0],[%i,5]),0);
  64 0$
  65
  66 (m : m . ctranspose(m),0);
  67 0$
  68
  69 (ok : true,0);
  70 0$
  71
  72 (for r in [complexfield, generalring, noncommutingring] do (
  73   ok : ok and checkcholesky(m,r)),ok);
  74 true$
  75
  76 /*-----------*/
  77
  78 (m : matrix([7-%i,0],[1-%i,5 + %i]),0);
  79 0$
  80
  81 (m : expand(m . ctranspose(m)),0);
  82 0$
  83
  84 (ok : true,0);
  85 0$
  86
  87 (for r in [complexfield, noncommutingring] do (
  88   ok : ok and checkcholesky(m,r)),ok);
  89 true$
  90
  91 /*------------*/
  92
  93 (m : matrix([matrix([1,0],[%i,10]), matrix([0,0],[0,0])], [matrix([1,2],[5,%i]),matrix([1,0],[2,8])]),0);
  94 0$
  95
  96 (matrix_element_transpose : 'transpose, matrix_element_mult : ".", 0);
  97 0$
  98
  99 (m : m . ctranspose(m),0);
 100 0$
 101
 102 checkcholesky(m, noncommutingring);
 103 true$
 104
 105 (matrix_element_transpose : false, matrix_element_mult : "*", 0);
 106 0$
 107
 108 /*---------------*/
 109 (ma[i,j] := kron_delta(i,j) + kron_delta(i,j+1)/2 + kron_delta(i,j-1)/2,0);
 110 0$
 111
 112 (m : genmatrix(ma,71,71),0);
 113 0$
 114
 115 checkcholesky(m,floatfield);
 116 true$
 117
 118 (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);
 119 0$
 120
 121 (matrix_element_transpose : 'transpose,0);
 122 0$
 123
 124 (matrix_element_mult : ".",0);
 125 0$
 126
 127 (m1 : cholesky(m.transpose(m), 'noncommutingring),0);
 128 0$
 129
 130 zeromatrixp(factor(m1.transpose(m1) - m.transpose(m)));
 131 true$
 132
 133 /*----clean up----*/
 134
 135 (remfunction(checkcholesky),0);
 136 0$
 137
 138 (matrix_element_transpose : false,0);
 139 0$
 140
 141 (matrix_element_mult : "*",0);
 142 0$
 143
 144 (remvalue(m,m1,ok),0);
 145 0$
 146
 147 (remarray(ma),0);
 148 0$
 149
 150
 151
 152
 153
 154
 155